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Abstract — Accurate channel models are of high importance 
for the design of upcoming mobile satellite systems. Nowadays 
most of the models for the land mobile satellite channel <LMSCt 
are based on Markov chains and rely on measurement data, 
rather than on pure theoretical considerations. A key problem 
lies in the determination of the model parameters out of the 
observed data. In this work we face the issue of state identification 
of the underlying Markov model whose model parameters are 
a priori unknown. This can be seen as a hidden Markov 
model jHMMt problem. For finding the maximum likelihood 
estimates of such model parameters the Baum-Welch 
<BW> algorithm is adapted to the context of channel modeling. 
Numerical results on test data sequences reveal the capabilities 
of the proposed algorithm. Results on real measurement data are 
finally presented. 



I. Introduction 

Satellite services to mobile users are experiencing a renewed 
interest thanks to the licenses granted for S-band usage for 
broadcast and interactive services |[T|-||3|. The underlying 
communication channel, referred to as land mobile satellite 
channel (ILMSCt . is characterized by strong variations of 
the received signal power. Obstacles in the propagation path 
between the satellite and the mobile terminal, such as buildings 
or trees may cause shadowing or even a complete blockage of 
the signal. With increasing frequency and decreasing elevation 
angle such events become more and more likely and strongly 
impact service availability. A further source of fading is due 
to multipath propagation: objects in the vicinity of the receiver 
are source of reflections that cause constructive or destructive 
interference. In the past several authors proposed Markov 
chain models to describe the behavior of the ILMSCI lU- 
JS). The modeling approach can be divided into two stages. 
First a Markov chain is set up to model slow transitions 
between different signal levels due to macroscopic effects such 
as blockage, shadowing, etc. In practice models with two or 
three states are common, but also a larger number of states 
is possible. Second, fast signal variations within each state 
due to multipath are taken into account assuming that the 
signal amplitude follows some specific distribution. To give 
an example, a Ricean distribution may be used to describe the 
signal amplitude in line of sight ( ILoSI l conditions, whereas the 
amplitude in a blockage state could be assumed to be Rayleigh 
distributed. 

Knowing the underlying channel model, a major issue 
consists in how to determine the model parameters out of 
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a sequence of measurement data. In literature there exist 
several approaches, most of them being rather simple and 
empirical. In lO the authors propose first to associate with 
each measurement sample a state of the underlying Markov 
chain. This association is done manually. Then, for each state 
the distribution of the associated samples is approximated 
by some known distribution by means of curve fitting. In 
m the weighted sum of some known distributions is fitted 
to the probability density function (p.d.f. i of the measured 



data. This gives the parameters for the distributions in the 
different states. Then, each sample is associated with a state 
by placing thresholds on the signal level. The thresholds are 
put according to the state probabilities from the fitting step. For 
highly overlapping distributions, this only works with limited 
accuracy, as we will show later. A more rigorous attempt is 
the technique in Q based on reversible jump Monte Carlo 
computation |[9l- It suggests fully blind estimation, making no 
prior assumptions on the number of states, nor on the specific 
distributions, allowing huge flexibility. However, this has the 
price of a significant increase in complexity and the resulting 
states and distributions often lack sufficient explanations in 
terms of underlying physical effects. 

Within this work we propose a further way to estimate 
the model parameters. It exploits the fact that the state iden- 
tification can be seen as a hidden Markov model ( IHMMI ) 
problem: out of the channel observation we would like to 
draw conclusions about the underlying Markov process that 
is not directly observable. A solution to this problem is given 
by the Baum-Welch dBWI ) algorithm, that has been widely 
used in other fields, such as speech or pattern recognition. 
An application to models of digital channels has already been 
provided in |10||. In the sequel our focus is on the ILMSCl 
We impose some constraints on the IBWI algorithm in order 
to improve its convergence and for sake of simplification. 
In particular it is well-known that its convergence properties 
depend on the initial model assumptions. Hence, unlike in Q 
we assume prior knowledge on the type of distributions and the 
distribution parameters (to be obtained by a preceding curve 
fitting step). Also, we fix the number of states in advance. 

The remaining part of the paper is organized as follows. In 
Section we recap the IBWI algorithm. Further we introduce 
a log-domain computation of the forward-backward metric of 
the IBWI algorithm and discuss some adaptations. Section |lll] 
reports the performance of the algorithm on test data, as well 
as on sequences of real measurement data. A comparison with 
the method in m is provided. 



II. Overview of the IBWI Algorithm 

Following the footsteps of ifTTl . we consider next the 
problem of associating an observed sample with a state of our 
IHMMI The IBWI algorithm can be applied to maximize the 
probability of a state given the entire observation sequence. 
Let us denote as Xt the state of the IHMMI at time t, and as 
r = (ri, r2, . . . , r„) the vector of n observations. The problem 
can be formalized as follows: given the vector of r, we are 
interested in locally calculating the probability of being in state 
i at time t, i.e. Pi{Xt — i\r}. We define 

fn.Xt (r, Xt = i) 



gtii) ^ PT{Xt ^ i\r} 



Mir) 
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For now, we focus on the joint p.d.f. in the enumerator. 
Under the IHMMl assumption, after dropping the subscripts for 
simplicity, this p.d.f.| can be rewritten as 



fir,Xt = z) = firlXt = z) ■ /(r^Vil^t = ^ 



at{i) bt{i) 

Here we used the shorthand r™ to denote the elements 
(rfc, rfe+i, . . . , Tuj) of the observation sequence r, with w > k. 
Further, referring to (|2|i, we define a forward metric at{i) and 
a backward metric bt{i)- It follows that 

at{i)bt{i) 



(2) 



9tii) 
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where the normaUzation by J2"Li'^t{i)bt{i) corresponds to 
fn.{r) in ([T]i and m denotes the number of states. Being pij 
the transition probability from state i to j, pi the probability 
of state i and fi{r) the probability density function given 
state i, the forward and the backward metric can be computed 
iteratively as 



at{i) 



l(.7>. 



bt+i{j)Ptjf]irt+i), 



(4) 
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with the initial metrics ai{i) = Pifi{ri) and 6„(i) = IVi. 

Figure [T] shows an excerpt of a state diagram for a Markov 
chain in the interval [t — l,t + 1]. The nodes of the trellis 
diagram at each time instance denote one of the three possible 
states, whereas the lines denote all possible transitions. Con- 
sider for example state 1. As indicated by the solid arrows, the 
metrics from all the nodes at i — 1, as well as t + 1 contribute 
to the calculation of the probability of state 1 at time t. The 
most likely state sequence can be determined by choosing the 
state with the highest probability at each time instance. 

Further, we may wish to estimate the probability of having 
a transition from state i at time t to state j at time t + l, given 
the observation r. This can be expressed as 



and it turns out that 
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Fig. 1. Calculation of the forward-backward metric for 3 states (linear 
domain). 



Even if the |BW| algorithm is in principle more general, we 
restrict ourselves to the simple case where the density func- 
tions fi{r), i = 1, . . . , TO, are perfectly known, whereas we do 
not have any knowledge about the transition probabilities pij, 
which we want to estimate. To do so, we chose some initial 
values for Q run the forward-backward algorithm and re- 
estimate the transition probabilities and the initial state 
probabilities pi according to the re-estimation formulae 
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Former values of pij and pi are replaced by the new estimates 
and the forward-backward algorithm is run again, leading to 
updated estimates, which are then fed-back. This process is 
iterated several times. The state identification step corresponds 
to the E-step of the expectation-maximization (lEMI i algorithm, 
where the model parameters are assumed to be fixed. The re- 
estimation step corresponds to the M-step, where the most 
likely model parameters are determined given the hidden state 
sequence. 

A final remark is related to the convergence of the algorithm. 
It is well known that the lEMI algorithm, as well as its special 
instance, the |BW| algorithm, increases the likelihood of the 
model iteration by iteration till it converges to a maximum 
value |[T2l . However, the algorithm may converge to a local 
maximum of the likelihood function, rather than to a global 
one. The convergence of the algorithm can be facilitated 
by limiting the set of a priori unknown model parameters. 
Alternatively, a set of various starting points can be considered. 



A. Log-domain Implementation of the \BW\ Alsoritlim 

Already for short observation sequences {n > 100) the 
forward-backward metric may get numerically unstable. As 

'in principle the choice is arbitrary. Nevertheless initial values not too far 
from the real values facilitate the convergence of the |BW| algorithm. Good 
starting points can be found in literature (e.g. in (4j)- 



a solution, for each t a normalization of the metric is usually 
performed iHl. Alternatively, a log-domain representation of 
the corresponding equations is proposed here. Let us define 
log-probabilities as 7t(i) \ngt{i), at{i) ^ Inat(i) and 
/3t(i) ^ ln6t(i), with In(-) being the natural logarithm. Then, 
(|3]l can be rewritten as 

7,(t) = a,(t) + - lnJ2T=i exp(a,(0 + 

V- J 

exp(«;i) 

Note that the last term can be solved in the log-domain 
by applying recursively the so-called max* operator (also 
known as Jacobi logarithm) that is defined as max*(Ki, K2) = 
ln(exp(ACi) +exp(K2)). Exploiting the identity 

max*(Ki, K2) = max(Ki, ki) + In (1 + cxp(— |ki — K2I)) 

and noticing that max*(Ki, K2, K3) can be recursively calcu- 
lated as max*(Ki, max*(K2, 'ta)), we have 

7,(t) = + - max*(a,(t) + 

i— l:m 

In a similar manner, using the shorthand 4>i(rt} ~ \nfi{rt), 
TTij — \npij and tt^ = Inpi we have that the recursions 



ai{t) = 0i(rf) + max*(aj(i- 1) +7rji) 

j = l:Tn 



with ai{l) 



i{ri) and 



= max + 1) + 7r,y + (t^jin+i)) 

j = l:m 



with /3i(n) = 0, Mi. Finally, for the re-estimation of the IBWI 
metrics we define Ctihj) ~ lnzt(i,j). Taking ^ we have 

Ct{i,j) = at{i)+TTij+(f>j{rt+i)+f3t+i{j)- 



max" I max * (at(i) + tt^ + ^j(rt+i) + A+i(j)) 

?'— l:m \ J — l:7n 

The estimation of the parameters proceeds as 



mux *Ct{i,j) 

t—l:n — l 



max 

t=l:ri-l 



max*Ct(«,j) 

J — 1 : rn 



max*Ci(i,j). 



while 



B. Restrictions on the \BW\ Alsorithm 

For modeling the lLMSCl applving the IBWI algorithm some a 
priori restrictions on the channel parameters have been applied. 
This is mainly motivated by two reasons. First, if reasonably 
good estimates of some channel parameters are available, their 
use may facilitate the convergence of the algorithm. Second, 
we consider important that the obtained results have a clear 
physical interpretation. To give an example, we would like 
states to be associated with different physical events, such as 
blockage of the signal or direct ILoSI During this work the 
following restrictions have been applied: 

• The type of distributions to be used has been fixed in 
advance. The original IBWI algorithm allows estimating 
the densities fi{r) iteratively as a mixture of Gaussian 
distributions |fTT|. It is however well-established that typi- 
cal propagation conditions (blockage or lLoSI for instance) 
can be accurately modeled by known distributions. 
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Fig. 2. Re-estimated transition probability pi2 vs. Bhattacharyya distance 
for the IBWI algorithm and threshold methods. 



Estimates of the distribution parameters are provided to 
the IBWI algorithm. Such estimates can be obtained for 
instance by a curve fitting step and are kept fixed through 
the IBWI re-estimation. Alternatively at each iteration the 
estimates could be refined, given the intermediate results. 
The number of states is fixed in advance corresponding 
to some physical events, such as total blockage of the 
signal by obstacles or iLoSI . 



III. Applications of the IBWI Algorithm 

The capabilities of the IBWI algorithm on different data 
sets are evaluated next. First we generate artificially a test 
sequence of samples and run iterative re-estimation. Knowing 
the original model parameters, our goal is to assess the 
quality of the re-estimations provided by the IBWI algorithm. 
A comparison with the commonly used threshold method and 
some derivatives is done. Second, the IB Wl algorithm is applied 
to data obtained from a measurement campaign. 

A. Application on Test Data Sequences 

Given a Markov chain with transition probabilities pij we 
generate a sequence of states x = {xi,X2, ■ ■ ■ , x„). For each 
state, an observation sample according to the associated |p.d.f.| 
is produced. For simplicity, we fix the number of states m to 
2. For state 1, we choose a Gaussian distribution with standard 
deviation cri = 0.2. To perform different tests, the mean value 
/ii ranges from 0.4 to 0.9. The Gaussian distribution associated 
with state 2 has mean /i2 = 1 and variance (T2 = 0.2. Since 
typically the ILMSCI is highly correlated @, we choose the 
state transition probabilities of the Markov chain 
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with corresponding state probabilities pi — 0.333 and p2 ~ 
0.667. The length of the state sequence (observation sequence) 
was set to n = 100000. 

Given the observations r = (ri, r2, . . . , r„) and the knowl- 



edge on the p.d.f. j, our iterative re-estimation algorithm is ran 



to determine the state sequence Xt, for i = 1 . . . n, as well as 
the transition probabilities pij and the state probabilities pi. 
It should be obvious that the closer the mean values of both 
Gaussian distributions are, the bigger shall be the deviation 
between the re-estimated state sequence x = (ii, X2, . . . , Xn), 
the associated re-estimated transition probabilities pij, as well 
as state probabilities pi and the actual values. To measure the 
distance of the two distributions /i (r) , /2 (r) associated with 
the two states, we make use of the Bhattacharyya distance 



B{h{r),h{r)) = - In / ^ h{r) ■ h{r) dr. 



For sake of comparison we also apply the threshold method 
to separate the states I?]- Samples below the threshold r are 
associated with one state, the ones above with the other We 
select the threshold r, such that the average error probability 



/i(r) dr +P2 f2{r) dr 
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is minimized. In addition, we assume a priori knowledge of 
the state probabilities (which could be e.g. provided by a 
previous curve-fitting step). For two Gaussian distributions 
with variances cti = (72 = cr, this yields 

.2 1 



("1 + A*2 



cr-' In £1 

P2 



2 (/i2 - Ml) 

Further, to suppress frequent state transitions (crossing of 
the threshold) we apply moving average filtering on the 
observation sequence. The span of the moving average is set 
to 10 or 20 samples. 

Figure |2] illustrates the estimated transition probability pi2 
versus the Bhattacharyya distance for the IBWI algorithm and 
the threshold methods with and w/o filtering. Despite close 
mean values of both distributions, the IBWI algorithm provides 
always accurate estimates for the state transition probabil- 
ity Pi2, whereas the threshold methods typically fail when 
B{fi{r), f2{r)) < 0.4. Empirically we found that best results 
for the threshold methods can be obtained with an averaging 
window span of 10 samples. It shall be noted however that in 
this case the resulting state probabilities deviate remarkably 
as illustrated in Table |I] It turns out that with increasing 
averaging window size even for Z5(/i(r), /2(r)) = 0.78 the 
estimated state probability pi is too low. Figure [3] depicts 
the share of wrongly labeled states in the estimated state 
sequence x, obtained through a comparison of the original 
state sequence x with x. Again the IBWI algorithm provides by 
far the best results, followed by the threshold methods with 
filtering. Note that for low Bhattacharyya distances the share 
of errors converges to 0.33 which corresponds to pi. 

B. Application on Measurement Data 

In fall 2008 a vast measurement campaign was carried out 
along the US East Coast in the framework of the European 
Space Agency (lESAb funded MiLADY project [13]. During 
the field trials the signal levels of the four satellite digital audio 
radio service (ISDARSt satellites were recorded with a mobile 
vehicular receiver A statistical channel model was derived 
out of the collected measurement data employing the IBWI 



TABLE I 

Pi FOR lBWi AND VARIOUS THRESHOLD METHODS: UNFILTERED (Tl), 
FILTERED WITH WINDOW SIZE 10 (TIO) AND 20 (T20) SAMPLES. 
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Fig. 3. Share of wrongly identified states vs. Bhattacharyya distance for the 
IBWI algorithm and threshold methods. A state at time t is considered to be 
wrongly identified if xt xt. 



algorithm. The proceedings are as follows: we first perform 
a curve-fitting step on the overall p.d.f. of r similar to lH. We 
obtain parameters of the distributions in the different states 
which serve as input for the IBWI algorithm. The resulting 
state probabilities are used to initialize ai{t). The curve 
fitting is performed using simulated annealing dSAI l fT4l . a 
fast meta-heuristic method for global optimization. In case 
the function to be optimized has several local maxima ISAI 
may overcome these and converge to the global minimum. 
Following literature, we chose three simple distributions to 
characterize the fast signal variations in the different states. 
The signal amplitude is assumed to follow a Rice distribution 
in case of direct iLoSl to the satellite. We associate a lognormal 
distribution with the shadowing state and assume that the 
signal amplitude in the blockage state is Rayleigh distributed. 
As second step, a preprocessing stage is required. The fast 
signal variations within a state are known to be correlated 
(see e.g. ifTSl ). However, (|2|i implicitly assumes independency 
among observation samples given a certain state. To comply 
with the independence assumption, the measurement data is 
down-sampled, taking into account the coherence time of the 
procesfl This leads the final observation r. Finally, given 
the three distributions and the observed sequence, the IBWI 
algorithm is applied as described in Section |ll] 

Let's consider a typical US urban environment with a 



^The spatial separation between samples after down-sampling was chosen 
to be 1 m in accordance with (4). 
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Fig. 4. Cui've fit on measurement data for urban environment and a satellite 
elevation of 30°. 



TABLE II 

Mean state duration D in meters and state probability pi for 

THE 3 PROPAGATION states. 
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satellite elevation of 30°. The solid line in Figure |4] shows 
the |p.d.f.| of the measured signal envelope, whereas the dashed 
lines with markers give the results of curve fitting using the 



three p.d.f. ; specified previously. The weighted sum of the 
Rice, lognormal and Rayleigh distributions is also plotted 



(dashed with diamonds) and turns out to be close to the p.d.f 
of the measured data. The Bhattacharyya distance between 
the lognormal (Rayleigh) and Rice (lognormal) distribution 
is 0.5 (1.1), thus posing challenges for state identification. 
Table HH gives the mean state durations Di = 1/(1 —pu) and 
state probabilities pi obtained with different state identification 
methods with minimum state duration set to 1 to. Results for 
the IBWI algorithm, the threshold method from ID, as well 
as the modified threshold method with a filter length of 10 
samples are shown. It can be observed that the IBWI algorithm 
and the threshold method yield the same state probabilities pi 
as obtained by means of curve fitting. However, the mean state 
durations obtained by the threshold method are very short. As 
illustrated in Figure |2] the threshold method tends to over- 
dimension pij, thus to under-dimension Di. If a prior filtering 
step is applied, the state durations become longer than the ones 
obtained with IBWI This is caused by an under-dimensioning 
of Pij for Bhattacharyya distances greater than 0.5 (c.f. Figure 
|2j. Here, the state probabilities are no longer preserved. 

IV. Conclusion 



of measurement data. The IBWI algorithm, allows estimating 
iteratively the hidden state sequence and the transition prob- 
abilities of the underlying IHMMI even for highly overlapping 
states. Especially in environments with frequent shadowing 
events conventional methods, such as the threshold methods 
may lead to inaccurate results on the state transition matrix 
(ISTMb of the hidden Markov process. Adaptations of the IBWI 
algorithm presented here guarantee numerical stability, as well 
as proper convergence at manageable complexity. Adaptations 
to channels different from the ILMSCl are possible and may be 
a matter of future investigation. 
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